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Abstract 



o 

^ In this work, the tridiagonal method lIU is used to distinguish between edges modes and area modes 



I to study the edge sites properties effect on edge localized states of semi-infinite zigzag 2D honeycomb 

^ graphene sheet. The results show a realistic behavior for the dependance of edge localized states of zigzag 

^ graphene on the edge sites properties which explaining the experimental results of measured local density 



of states at the edge of graphene [2], while at the same time removing the inconsistence between the semi- 
conductor behavior found in the experimental data for fabricated GNRs [31 |4l and the expected theoretical 



semi-metallic behavior calculated without considering the edge properties effect on the edge localized states 
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I. INTRODUCTION 



It is known that the 2D materials with zigzag edged nanoribbons of 2D honeycomb lattice 
structure has peculiar flat localized edge states at the Fermi level [HllSHHll, which is a result of the 
zigzag geometry effect on the particles hopping flow in its edges sites. 

These edge states are known to be important due to their efi'ect on the electronic properties and 
as a consequence in a variety of future applications of the famous Zigzag Graphene Nanoribbons 
(ZGNR) yj|8|. The edge localized states are, in general, dependent upon the ribbon size and purity 
of the sample [[BISll. From a theoretical point of view, the edges states depends on the probability 
for an electron to hop from a site in the edge to a bulk site, or to impurity site in the neighborhood. 
The edge atoms have a different coordination number from the bulk atoms, this leads to a diff"erent 
hopping parameter between the edge atoms and the bulk one. Such a diff"erence is not usually 
considered in previous calculations for that ZGNR edge localized states with different approaches 
[ISHHl, these calculations show inconsistencies with experimental results for all fabricated GNRs 
that has semiconductor behavior [OIH which is a consequence of the absence of the flat edge states 
at Fermi level. 

In [HOl we found that the tridiagonal method has an advantage in studying the edge properties 
effects on the edge localized states due to its ability to separate the edges modes from area modes 
in case of the 2D square lattice. 

Therefore in this work, the tridiagonal method is used to study the effect of edge sites properties 
on the edge localized states of semi-infinite zigzag 2D honeycomb sheet as study case. The similar- 
ity between semi-infinite ZGNR and semi-infinite antiferromagnetic as both two-sublattice struc- 
ture guide us to follow the tridiagonal method steps used in study the surface modes of Heisenberg 
antiferromagnetic [fm[T2| . The method allow us also to study the effect of impurities introduced 
substitutionally in impurities localized states of the semi infinite ZGNR. 

II. THEORY FOR EDGE STATES AND IMPURITY STATES 

The structure of semi-infinite ZGNR is a honeycomb lattice of carbon atoms with two sublat- 
tices denoted as A and B. The geometry of a graphene nanoribbon with zigzag edges is shown in 
Figure [T| where the system is infinite in the x direction and has 2N rows of carbon atoms in the 
3; direction. To be considered a ribbon, A^^ is a finite integer but here we will extend the study for 
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FIG. 1. Geometry of a graphene ribbon with zigzag edges. The black (gray) dots are the sublattice A (B) 
atoms, where A{B) sublattice type are labeled by index n{n') {= 1, 2, • • • ,N) and the white dots show a row 
of impurities. Figure taken from [1131 . 

the the semi-infinite case where N ^ oo. The A{B) sublattice type lines are labeled with index 

n{n') where n(n') = 1, 2, 3, The impurities (which may be silicon or boron, for example) are 

introduced substitutionally along two different rows of atoms parallel to the x axis. The impurity 
lines, which preserve the translational symmetry in the x direction, may be any distance apart in 
the ribbon. 

TABLE I. Nearest neighbor hopping matrix elements for the zigzag graphene nanoribbon 

Parameter Zigzag 
jS 2fcos( V3^;ta/2) 

r t 



Following a microscopic approach in terms of a tight binding model Hamiltonian [8J with 
neglecting the next nearest neighbor and with follow some recent work for impurities in graphene 
ribbons [fT^ the Hamiltonian becomes 

H = -Y,tiM]bj + aib]) (1) 

where a] (or a,) creates (or annihilates) an electron on the sublattice A site, and Z?j(or bj) does the 
same for the sublattice B site, while tij is the nearest-neighbor hopping energy between sublattices. 
In the pure material the hopping energy is denoted by t and its value is known JSl to be ~ 2.8 eV. 
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Taking into account the translational symmetry in x direction, a Fourier transform is made to 
rewrite Equation ([T]) in a wavenumber representation q^c in the x direction, and the rows are labeled 
n and n'. The Hamiltonian becomes 



qx,nn' 



The hopping amplitude factors Tnn'iqx) for the zigzag structure have the form 



(2) 



'^nn'^Qx) ^ 



2 COS 



rv3 



-q^a 



(3) 



or 



(4) 



where the assignment of upper or lower signs depends on the sublattice type sequence for rows n 
and n' (see Appendix [A]). The definition of /3 and y is given Table |lj Now we use the equation of 
motion iMX/dt = [X, H] for any operator X for the creation and annihilation operators of each 
row. Taking fi = I and assuming that the modes have a time dependence like exp[-ia)(qx)t], we 
obtain 2N coupled equations: 



<^{qx)bq^,n' = ^Tn'n{qx)aq,,n'- 



(5) 



Expanding and rearranging the Equations ([5]) (see Appendix |A]) such that coupled equations 
between sublattice A and sublattice B operators, could be written in the following form 

[ojHqx)-[0'+y'')] 



Equations ([6]) could be written in the following supermatrix equation [fTTTl 



(6) 



'' An + AAn On ^ 
Bn In 



0, 



(7) 
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where is the null matrix, /jv the identity matrix, aNib^) operator column vector, and 

^^-10 
-1 ^ -1 
An= 0-1^-10 
0-1^-1 



and 



B 



N 



77 

/I 77 

/I 77 

A J] 



(8) 



(9) 



the elements of above matrices are defined by 



T] = 



-7 



A = 



(10) 



fir io(,qx) 
The edge properties have been separated from the area "bulk" properties of ZGNR by forming 

the matrix IsA^. To simplify the calculations we consider only putting one or two impurities lines 

in rows numbers uq and h'q of sublattice A such that their properties could be separated form area 

properties in the same way the edge properties separated before by including them in the matrix 

AA^. In this case the matrix AA^ has the following form 

A, A, 

A, 

0000 000 

A7„o 

Ano A,„„ 



AAn = 































A/„. 













'"0 





















A„' A/,,' 

"0 "'0 



A 



In'g 












(11) 
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the elements of AA/v matrix are defined by 



fir ' 



5"0 



;0y 



(12) 



/3y '"0 

where the edge hopping te, the first impurities line hoping and the second impurities line hoping 
t„'^ have replaced the ZGNR interior area sites hoping t in the definition of /3 and y in Table |l] to 
obtain the edge and impurities counterpart. 

Follow the steps of Heisenberg antiferromagnetic case [fTTll and the algebra of block matrices 
[fT4|. one define the supermatrix G 



G = 



^ An On^ 
\ Bn In j 



' (A^)-i On ^ 
, -Bn^An) ' In 



(13) 



Multiplying Equation ( [13] ) in EquationQ, we get the following 

{In + {AnY^ ^AN)aN = 
(-BMNy^^N)aN + bN = 

Define the following matrix 

Dn = In + (AnT^AAn 
The matrix Dn could be written in the following partition form (see Appendix |B]) 

= 



(14) 



(15) 



(16) 



where O is a square null matrix, / a square identity matrix, S a square submatrix of Dn, and Q is 
square submatrix of Dn with dimension of n[, + 1 x + 1. 

The elements for the inverse of tridgional matrix An [[TTl[T5] - [T8l . i.e. the matrix (An)'^ is given 
as following: 

yi'i+'n y\n-~m\ 

i(ANr\m = (17) 

X + X ^ 

where jc is a complex variable such that \x\ < 1 and x + x"^ = ^. As mentioned in [[TOl . the 
values of x should satisfy the following boundary and physical conditions [fT9ll . The area modes 



are oscillating waves inside the nanoribbon, which requires that x must be imaginary exponential 
X = e'*"''^ with |;c| = 1. From the definition of ^ and x parameters, the dispersion relation for the 
area band is given by 



^ = x + x-^ = + = lcos{qyall) 



(18) 



coB(qx,qy) = ±hU +4cos2 



+ 4 COS 



cos 



This expression for the 2D area band for the zigzag nanoribbons is very similar to the extended 
graphene electronic dispersion relation given in dH. This expression also shows the same general 
features of graphene band structure ||8l. 

While the edge modes are localized on the edge and they are decaying exponentially inside the 
nanoribbon, which requires that x must be real and less than 1 for edge modes. The edge modes are 
obtained by requiring the determinant of the coefficients for operator column vector to vanish 
|III1[I91|201: 



\Dn\ = det 



Q 





s 


I 



(19) 



using the rules for obtaining the determinant of partitioned matrices [[T4|. Equation ([19]) become 

\Dn\ = \Q\\i-sq-'o\ = \q\, 



(20) 



the localized edge and impurities states, i.e. the edge and impurities dispersion relations, for the 
semi-infinite zigzag are obtained by taking the limit of Equation pO] ) as A/^ ^ oo 



lim ID^I = lim 121 = \Q\. 

N—nx N—>co 



(21) 



In the case of taking only the edge properties effect on the edge localized states, i.e there is no 



0, 



any impurities lines inside the sheet. Equation ( [2T| ) become 

det(e) 
which give 

A^x^ - 2 A,/ - (2A^ + A,)x^ - (1 - 2A,)x^ + (A^ + A,)x -1=0. 



Aj/A, +Aj]A, + 1 Aj/A, 

A22^Ae + A22AS Aj^Ag + 1 



(22) 



(23) 
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In case the interaction of the edge with the interior sites is not alfected with the edge sites 
properties, i.e. = 0, Equation pT] ) become 



det(A) = 



A\lK + 1 
A-|A, 1 



= A1\K + \=Q (24) 



which is the same expression obtained for Heisenberg antiferromagnetic [[TTI . 
III. RESULTS 

Figures [2] and |3] show the edge localized states of semi-infinite zigzag graphene sheet for dif- 
ferent edge hopping to area hopping ratios calculated with A, = using Equation [24] for Figure |2] 
and calculated with A^. using Equation [22] for Figure [3] The Figures show that the dispersion 
of the edge localized states depends on both the edge sites hopping properties and their effect on 
the interaction with the interior sites in the zigzag sheet. 

The Figures begin with edge hopping equal to zero, which could be done by saturating the 
carbons atoms on the edge. In this hopping value both calculations (with and without A J result 
in an extended flat localized edge state through the whole Brilloin zone at Fermi level ojpit = 0. 
This is due to the localized edge wave functions which agree with density-functional theory (DFT) 
calculations for finite ribbon []2T]| . 

As the edge hopping increases from zero, the dispersion of the edge localized state begin to 
have a percentage of it laying at Fermi level tOf = 0, and the remaining percentage liftoff from the 
Fermi level. The percentage of the edge dispersion that laying at Fermi level is important to the 
electronic properties of the zigzag edged graphene nanoribbon. 

To study the variation of the edge dispersions percentage laying at Fermi level with the chang- 
ing of the edge hopping and in the same time to easily clarify the data displayed in Figures ]2] and 
[3] we define the following three parameters for each given edge hopping value. The first one is 
the Relative Localized Density of States near Fermi level (FRLDOS), which is calculated compu- 
tationally by counting the total number of points in the localized edge dispersion with co/t < 0.2. 
The second parameter is the Total Relative Localized Density of States (TRLDOS), which is cal- 
culated computationally by counting the total number of points in the localized edge dispersion, 
which is a relative measure of its total density of states. The third parameter is the difference 
between FRLDOS and TRLDOS which represent the liftoff percentage of the edge localized state 
dispersion from Fermi level (LOFRLDOS). 
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FIG. 3. Edge localized states for edge with differen^edge hopping calculated with + 0, the shaded band 
represent area modes continuum. 




FIG. 4. The variation of TRLDOS, FRLDOS, and LOFRLDOS as a function of edge hopping form to 
2.9. (a) for A, = and (b) for A, t 0. 

Figure|4] shows the variation of TRLDOS, FRLDOS, and LOFRLDOS for edge localized states 
as a function of edge to area hopping ratio form to 2.9 with increment of 0.1, the (a) sub Figure 
represent the calculation with = 0, while (b) sub Figure represent the calculation with 4^ 0. 

With the help of the comparison between Figure^a) and Figure[2]and the comparison between 
Figure Qb) and Figure [3j the variation of the edge localized states dispersion with different edge 
hopping values is described as follow: 

First, for the dispersion of localized edge state calculated using A^ = 0, starting with edge 
hopping equal to zero the Figure |4]^a) shows that FRLDOS is equal to TRLDOS which means 
that the edge localized state lays completely in the Fermi level while LOFRLDOS is zero. In 
Figure |2] it is shown as a flat localized edge state at Fermi level cop/t = extended through the 
whole Brilloin zone as described above. For 2 increment in the edge hopping there is no change 
in the value of FRLDOS and it is still equal to TRLDOS with LOFRLDOS is equal to zero, which 
mean that dispersion did not change from as zero hopping. Beginning from edge hopping equal 
to 0.3 the value of FRLDOS begin to decrease and LOFRLDOS begin to increase while TRLDOS 
keep constant which means that some of dispersion left off from Fermi level as shown in Figure 
[2] as a left off near qxajn = +0.5. As edge hopping reach 0.5 the TRLDOS still keep constant 
while FRLDOS still decreasing and LOFRLDOS still increases which means that more dispersion 
left off from Fermi level as shown in Figure |2] as increase in the left off near qxa/n = +0.5 and 
around qxa/n = +0.0. At edge hopping 0.8 the TRLDOS drop quickly to small value and it is 
equal to LOFRLDOS while FRLDOS become zero which means that most the edge localized 
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dispersion disappear which shown in edge hopping 0.9 in Figure [2] as small edge localized points 
near the intersection of the area band segments. When edge hopping is equal to interior hopping, 
the edge localized states completely disappear. Figure |4]^a) shows that FRLDOS is equal to zero 
in edge hopping range 0.9 - 2.9. As the edge hopping increases from 1 to 2.9, LOFRLDOS is 
equal to TRLDOS. The LOFRLDOS increases very quickly to large value with increasing of edge 
hopping from 1 to 1.3, and then keep constant in the edge hopping range 1.3 - 1.7 after that 
range LOFRLDOS decrease slowly. This behavior is shown in Figure [2] as increasing in the edge 
localized states above the area band. 

In the second case, the dispersion of localized edge state calculated using A, 4^ 0, starting 
with edge hopping equal to zero the Figure |4]^b) shows that FRLDOS is equal to TRLDOS which 
means that the edge localized state lays completely in the Fermi level while LOFRLDOS is zero, 
in Figure [3] it is shown as a flat localized edge state at Fermi level io^li = extended through the 
whole Brilloin zone as described above. At the edge hopping equal to 0.1 the values of FRLDOS 
and TRLDOS begin to decrease and LOFRLDOS begin to increase which means that some of 
dispersion left off from Fermi level. In the edge hopping range from 0.2 to 1.4, the FRLDOS 
become zero and TRLDOS become equal to LOFRLDOS. In edge hopping range 0.2 to 1, the 
LOFRLDOS decreases to a zero value, which is shown in Figure [3] as decreasing in the edge 
localized states dispersion, and at the same time the edge localized state shifting up in the energy. 
It is very important to note that at edge hopping 0.5, the edge localized states dispersion become 
very similar to the famous peculiar edge localized state for graphene zigzag nanoribbon jSHSl but 
here shifted from Fermi level due to the edge hopping properties. When edge hopping is equal to 
interior hopping, the edge localized states completely disappear. The LOFRLDOS increases very 
quickly to large value with increasing of edge hopping from 1 to 1.2, and then remains constant 
in the edge hopping range 1.2 - 1.4. After that range, the LOFRLDOS increase to a peak at edge 
hopping 1.6 and then it begins to slowly decrease until edge hopping 2.9. Starting from the edge 
hopping 1 .4, the LOFRLDOS begins to slowly increase with increasing the edge hopping until the 
edge hopping reach 2.3, then the LOFRLDOS begins converge to a nearly constant value. While 
the TRLDOS increase to a peak at edge hopping 1 .6 and it then begins to slowly decrease with 
increasing the edge hopping. The behavior of the three parameters is displayed in Figure [2] as a 
change in the edge localized states around the area band. 

The efi'ects of both the impurities hopping and the impurities line position on the impurities 
localized states have been calculated using Equation |2T] and the results are shown in Figures [5} 
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The impurties line position: 2 
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The impurties line position: 4 
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The impurties line position: 3 
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The impurties line position: 5 
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FIG. 5. The variation of TRLDOS, FRLDOS, and LOFRLDOS as a function of impurities hopping form 
to 2.9 for the impurities line in sublattice A at positions (a) 2, (b) 3, (c) 4, (d) 5. 



This calculation was done considering the edge hopping is equal to 1 . 

Figures [5] show the variation of TRLDOS, FRLDOS, and LOFRLDOS as a function of im- 
purities hopping form to 2.9 for the impurities line in sublattice A at positions (a) 2, (b) 3, (c) 
4, (d)5. It is clear that the dispersion of the impurities localized states for impurities line at sec- 
ond raw of the sublattice A is very similar to the dispersion of the edge localized states described 
above. Beginning from position three in the sublattice A, the dispersion of the impurities localized 
states becomes completely different and nearly independent on the impurities Une position in the 
sublattice A. 

Figures [5j a), (b), and (c) show that for the impurities hopping range form to 0.9 the FRLDOS 
is equal to TRLDOS which means that the impurities localized state lays completely in the Fermi 
level while LOFRLDOS is zero. In Figure [3] it is shown as a fiat localized edge state at Fermi 
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level coplt = extended through the whole Brilloin zone as described above. At the impurities 
hopping equal to 1 the FRLDOS switches to zero while LOFRLDOS changes to the value that 
keeps TRLDOS constant. This means that the impurities localized state is completely lay off 
from Fermi level with keeping its density of states constant. As impurities hopping from 1 to 
1.3 the LOFRLDOS rapidly increases to a peak in its value, it then decreases until the impurities 
hopping reach 1.7 then the LOFRLDOS begin converging to nearly a constant value. Likewise, as 
impurities hopping from 1 to 1.7 the FRLDOS increase until the impurities hopping reach 1.7 the 
FRLDOS begin converging to nearly a constant value lower that LOFRLDOS value. While the 
TRLDOS has fast increases in edge hopping range from 1 to 1.3 and then it begin converging to a 
nearly high constant value. 




q a/n 



FIG. 6. Edge and impurities localized states the black dots for edge and two lines of impurities at sublattice 
A, the shaded band represent area modes continuum. The Edge hopping is tg/t = 0.5, the first impurities 
line position is «o = 4 with impurities hopping = 1.5, and the second impurities line position isriQ-l 
with impurities hopping tnjt - 0.3. 
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IV. DISCUSSION AND CONCLUSIONS 



In this work the tridiagonal method was used to study the effect of edge and impurities sites 
properties on their locahzed states in semi-infinite zigzag edged 2D honeycomb sheet. It is found 
that the tridiagonal method calculations provide us with two possibility to study the effect of the 
edges sites properties on their localized states. In the first one, the interaction of the edge with the 
interior sites is not affected with the edge sites properties, i.e. = 0, and in the second one the 
interaction of the edge with the interior sites is affected by the edge sites properties, i.e. 4^ 0. 

The results of the case A^ = show that the edge locahzed states dispersion has qx dependance 
of the hopping in ID chain [|22ll in low edge hopping values, which reflect the completely isolation 
of edge hopping from interior one. This behavior is away from the expected behavior from tight 
binding model results and the explanation given in [HI. 

The results of the case A^ 4^ show that the edge locahzed states dispersion has dependance 
of the hopping in 2D honeycomb but shifted in the energy due to the edge hopping properties and 
especially at edge hopping 0.5, the edge localized states dispersion become very similar to the 
famous peculiar edge localized state for graphene zigzag nanoribbon [|5]-|8l but here shifted from 
Fermi level, which reflect the importance of edge sites properties on the edge hopping with interior 
sites. In this case the behavior is agree with the expected behavior from tight binding model results 
and the explanation given in [[1]. 

In case of the effects of both the impurities hopping and the impurities line position on the 
impurities localized states, the results show that second raw of the sublattice A has very similar 
edge localized states dispersion, and after that position the impurities localized states become 
nearly independent on the impurities line position in sublattice A. Also, the results show that 
introducing impurities in any position of the sheet will produce impurities localized states at Fermi 
level in all impurities hopping properties, which affecting the electronic properties of the sheet. 

The model could be used to study the effects of hopping properties interaction for the edge and 
the two separated lines of impurities at sublattice A in the edge and the impurities localized states 
as shown in Figure |6] 

The equivalent between the obtained mathematical expressions for edge states in case of zigzag 
graphene and that obtained for surface spin wave in case of Heisenberg antiferromagnetic [fTT]| 
reflect their equivalent from geometrical and topological point of view. In the same time it show 
that result is applicable to the magnetic and the 2D materials have the same geometrical and 
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topological structure. 

Finally, the results of considering the interaction of the edge with the interior sites is affected by 
the edge sites properties, i.e. A , 4^ 0, show a realistic behavior for the dependance of edge localized 
states of zigzag graphene on the edge sites properties which explaining the experimental results of 
measured local density of states at the edge of graphene [2], and in the same time removing the 
inconsistence between the semiconductor behavior found in the experimental data for fabricated 
GNRs [|3l in and the expected theoretical semi-metallic behavior calculated without considering 
the edge properties effect on the edge localized states [EHSl. 
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Appendix A: Rearranging Equations Q 



In this appendix we list the rearrange steps of Equations (|5]) to obtain Equations ([6)). 



i^{qx)bq^,n' = y^Tn'n{qx)aq„n'- 
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Expanding Equations (|5]) using the sublattice indexes n and n' 

Oj(qx)aqx,n = ybq,,n'-\ + /3bq^,n' 



oj(qx)bq,,i = /3aq^,i + yaq^2 

Oj{qx)bq,,2 = Pciq^,2 + yaq,,i, 

oj(qx)bq,,3 = Pag„3 + ycigxA 

<^{qx)bq„n' =Paq„n + ydq^n+l 



divide by a>(qx), and rearrange we get 



Oj(qx)bq^,n' = fiaq,,n + yaq^^n+l 



Oj(qx)bq„n'-l =l3aq^^n-\ + yClq^n 



Oj(qx)aq^,n = ybq,,n'-l + fibq.^n' 



, , Paq„n-\ + yaq„n ^ficiq^n + yClq^n+l 

(^(qx)aq,,„ = y -7^7^ 



OJ{.qx) 



ii>{qx) 



(qx)aqx,n = y{fiaq^^n-\ + Tfl^^n) + l^ifittq^^n + yaq„n+\) 
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lead to 



= 



Appendix B: The partition of matrix 



In this appendix we list the steps of partitioning the D^^ matrix. 



4-1 
^11 


4-1 
^12 


4-1 • 

^13 




Ino 


4-1 


^21 


4-1 
^22 


4-1 • 

^23 


•• 4-1 


4-1 


4-1 
^2no+l 


4-1 
^31 


4-1 
^32 


4-1 • 

^33 


^3no-l 


4-1 


4-1 



nA HrJ. n„: 



-1 




A~\ 


K- 


1 




-1 




A-\ 




1 




-1 






3«o- 


1 


3«o 



A-1, , A-i, A 
A / , , A / / A 



which is of dimension N xN. 



AA^ = 



A. 


A. 





























A. 













































































A/„o 





























A/no 


A„o 


A/no 





























A/„o 













































































A/„' 





























A/„' 

'"0 


A„' 

"0 


A/„' 
"'0 





























A/„' 










































which is of dimension N x N. 
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Dn = In + (AnT'AAn (B3) 
therefore the Dj^ dimension is N x N with the following elements: 



Ai 






= AJi^A, + 6i2 


ItlQ—V 


= y4- ^ Ar„„ + (5,„„ 1 
(no ^"0 "mo-i 














Din' 


= A:!_iA,„, + A:1a„. + A7„1,iA,„. + cJ,„. 


Din'g+1 






= if both j and / are not equal to either 



Dn 


D12 


Dn ■ 


• Dino-l 


Dlno 


Dlno+1 • 






D 


D21 


D22 


D23 • 


■ -02,70-1 


D2na 








D 


D31 




As • 


• D^no-l 


Dino 


Ano+1 ■ 






D 



Dn^jl DnQ2 Dn^'i ' ' ' Dn^^na-l -Onnm Dn 



'nofio-l ^nono ' ' ' Dn^n'^-l Dn^n'^ D„f^„'^+i 

Dn'gl Dn'^2 Dn'^3 ' ' ' Dn'^UQ-l Dn'^no ' ' ' '^"o^o"! '^"o''o ^"o^O"*"^ 

which give the following partition of D^ matrix 



D 



N 



Q 


o) 


s 


n 



(B4) 



where O is a square null matrix, / a square identity matrix, S a square submatrix of Df^, and Q is 
square submatrix of D^ with dimension of + 1 x nQ + 1. 
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